2-2

Hexagonal Coordinate System

Nicholas Wheeler
20 January 2012

Pre-introduction

While this work began as an attempt simply to discover an efficient way to coordinatize/identify/distinguish the nodes of a hexagonal grid, it proceeds to a sketch—made possible by that accomplishment—of the Markovian approach to random walks on such a grid. I am led to the conclusion that such an approach would require the resources of a supercomputer.

Introduction

Here is the a portion of the unbounded hexagonal lattice on which, in recent work, I have constructed weighted random walks:

The nodes of the fundamental cell stand at noon, 2:00, 4:00, 6:00, 8:00 & 10:00 o'clock on the unit circle (center at the origin):

P[n_]:={Re[Exp[I((π/2)-n(π/3))]],Im[Exp[I((π/2)-n(π/3))]]}

Graphics[{
{PointSize[0.03],Point[P[0]],Point[P[1]],Point[P[2]],Point[P[3]],Point[P[4]],Point[P[5]]}}, Axes->True, PlotRange->{{-1.1, 1.1},{-1.1, 1.1}}, Ticks->{{-((Sqrt[3])/2),((Sqrt[3])/2)}, {-1,-(1/2), (1/2),1}}]

My problem is to assign coordinates to nodes in such a way that (1) the coordinates of nearest-neighbor nodes are easily identified and (2) it becomes possible to associated nodes in a natural way with the elements of a Markov matrix and of the stochastic vectors on which those matrices act.

Basic idea: red/blue nodal lines and their points of intersection

We observe that each red (ditto blue) node lies at the intersection points of a pair of sloping node rows. We undertake, therefore, to describe those linear rows.

Generally, the line on the (x, y)-plane that passes through the points ((Subscript[a, 1]),(Subscript[a, 2])) and  ((Subscript[b, 1]),(Subscript[b, 2])) is defined implicitly by the requirement that

Det[(({{1, x, y}, {1, (Subscript[a, 1]), (Subscript[a, 2])}, {1, (Subscript[b, 1]), (Subscript[b, 2])}}))]

%==x((Subscript[a, 2])-(Subscript[b, 2]))-y((Subscript[a, 1])-(Subscript[b, 1]))+((Subscript[a, 1])(Subscript[b, 2])-(Subscript[a, 2])(Subscript[b, 1]))//Simplify

To obtain the lines that pass with positive/negative slope through the node at "noon" we construct

Det[(({{1, x, y}, {1, 0, 1}, {1, -((Sqrt[3])/2), -(1/2)}}))]
Det[(({{1, x, y}, {1, 0, 1}, {1, +((Sqrt[3])/2), -(1/2)}}))]

ContourPlot[{((Sqrt[3])/2)+((3 x)/2)-(((Sqrt[3]) y)/2)==0,-((Sqrt[3])/2)+((3 x)/2)+(((Sqrt[3]) y)/2)==0},{x,-2,2},{y,-2,2},Axes->True, AxesLabel->{x,y}]

We are interested in the x-translates of those lines, where the elemental translation carries

and has length (Sqrt[3]). We re led thus to define

Rup[x_, y_,n_]:=((Sqrt[3])/2)+((3 (x-n(Sqrt[3])))/2)-(((Sqrt[3]) y)/2)
Rdn[x_, y_,n_]:=-((Sqrt[3])/2)+((3 (x-n(Sqrt[3])))/2)+(((Sqrt[3]) y)/2)

where the function names suggest "Red up-slope" and "Red down-slope."

The red nodes of the fundamental hex-cell lie at three of the four points where these lines (two up-sloping, two down-sloping)

Join[Table[Rup[x,y,n]==0,{n,0,1}],Table[Rdn[x,y,m]==0,{m,-1,0}]]

intersect. We plot those lines

BasicRedLines=ContourPlot[{((Sqrt[3])/2)+((3 x)/2)-(((Sqrt[3]) y)/2)==0,((Sqrt[3])/2)+(3/2) (-(Sqrt[3])+x)-(((Sqrt[3]) y)/2)==0,-((Sqrt[3])/2)+(3/2) ((Sqrt[3])+x)+(((Sqrt[3]) y)/2)==0,-((Sqrt[3])/2)+((3 x)/2)+(((Sqrt[3]) y)/2)==0},{x,-3,3},{y,-2,2}, AspectRatio->Automatic,Axes->True]

The "noon," "4 o'clock" and "8 o'clock" nodes have coordinates given respectively by

{x,y}/.Solve[{Rup[x, y,0]==0,Rdn[x, y,0]==0},{x,y}][[1]]

{x,y}/.Solve[{Rup[x, y,1]==0,Rdn[x, y,0]==0},{x,y}][[1]]

{x,y}/.Solve[{Rup[x, y,0]==0,Rdn[x, y,-1]==0},{x,y}][[1]]

which we superimpose upon the preceding figure:

BasicRedNodes=Graphics[{Red,PointSize[0.03],Point[{{0,1},{((Sqrt[3])/2),-(1/2)},{-((Sqrt[3])/2),-(1/2)}}]}];

Show[{BasicRedLines, BasicRedNodes},Axes->True]

The down/up-sloping blue-node lines are got from the up/down-sloping red-node lines by this simple adjustment (reflection in the x-axis):

y⟶-y

Whence the constructions

Bdn[x_, y_,n_]:=((Sqrt[3])/2)+((3 (x-n(Sqrt[3])))/2)+(((Sqrt[3]) y)/2)
Bup[x_, y_,n_]:=-((Sqrt[3])/2)+((3 (x-n(Sqrt[3])))/2)-(((Sqrt[3]) y)/2)

BasicBlueLines=ContourPlot[{Bdn[x,y,0]==0,Bdn[x,y,1]==0,
Bup[x,y,0]==0,Bup[x,y,-1]==0},{x,-3,3},{y,-2,2}, AspectRatio->Automatic,Axes->True]

{x,y}/.Solve[{Bup[x, y,0]==0,Bdn[x, y,0]==0},{x,y}][[1]]

{x,y}/.Solve[{Bup[x, y,0]==0,Bdn[x, y,1]==0},{x,y}][[1]]

{x,y}/.Solve[{Bup[x, y,-1]==0,Bdn[x, y,0]==0},{x,y}][[1]]

BasicBlueNodes=Graphics[{Blue,PointSize[0.03],Point[{{0,-1},{((Sqrt[3])/2),(1/2)},{-((Sqrt[3])/2),(1/2)}}]}];

Show[{BasicRedLines, BasicBlueLines, BasicRedNodes, BasicBlueNodes}]

We have actually no interest in the red/blue nodal lines themselves (which were introduced only as a coordinating device), but only in their points of intersection. Accordingly,  we use

{x,y}/.Solve[{Rup[x, y,m]==0,Rdn[x, y,-n]==0},{x,y}][[1]]

to construct an {m.n}-indexed array of points (red nodes)

PR[m_,n_]:={((Sqrt[3])/2) ( m- n),1-(3/2) ( m+ n)}

Simplify[PR[m,n]=={((Sqrt[3])/2) ( m- n),1-(3/2) ( m+ n)}]

of which we plot a sample:

RedPopulation=Flatten[Table[PR[m,n],{m,-5,5},{n,-5,5}],1];

RedPopFig=Graphics[{Red, PointSize[0.02],Point[RedPopulation]}, Axes->True, AxesLabel->{x,y}, Ticks->False]

Similarly,

{x,y}/.Solve[{Bup[x, y,m]==0,Bdn[x, y,-n]==0},{x,y}][[1]]

PB[m_,n_]:={((Sqrt[3])/2) ( m- n),-1-(3/2) ( m+ n)}

Simplify[PB[m,n]=={(1/2) ((Sqrt[3]) m-(Sqrt[3]) n),-((2 (Sqrt[3])+3 (Sqrt[3]) m+3 (Sqrt[3]) n)/(2 (Sqrt[3])))}]

BluePopulation=Flatten[Table[PB[m,n],{m,-5,5},{n,-5,5}],1];

BluePopFig=Graphics[{Blue, PointSize[0.02],Point[BluePopulation]}, Axes->True, AxesLabel->{x,y}, Ticks->False]

Show[{RedPopFig, BluePopFig}]

In construcing the red (ditto blue) nodes in this display we let m and n range independently on

so there are (11^2)= 121 red nodes, an equal number of blue nodes, and 242 nodes in total. 

If those nodes were connected—each to its to nearest neighbors to construct complete or (at the edges) partial hexagons—the top red and bottom blue nodes would become disconnected, leaving a connected graph with 240 nodes. Markov matrices that act on 240-dimensional stochastic vectors have (240^2)= 57,600 elements. Since each node is connected to (at most) three others, such matrices would be "sparce," so computational management may not be out of the question.

Nearest neighbors

Which three of the blue nodes PB[p, q] are nearest neighbors  of the red node PR[m, n]? 
Which three of the red nodes PR[p, q] are nearest neighbors  of the blue node PB[m, n]? 

In previous work I gave names as indicated in the following figure

to the displacement vectors

R1={Cos[(π/2)],Sin[(π/2)]};
R2={Cos[(π/2)-4((2π)/12)],Sin[(π/2)-4((2π)/12)]};
R3={Cos[(π/2)-8((2π)/12)],Sin[(π/2)-8((2π)/12)]};

that proceed ● ⟶ ●, and

B1=-R1;
B2=-R2;
B3=-R3;

to  those that proceed ● ⟶ ● :

Specifically,

R1
R2
R3

and

B1
B2
B3

A little experimentation  establishes that those vectors can be described

R1==Simplify[PB[m-1,n-1]-PR[m,n]]
R2==Simplify[PB[m ,n-1]-PR[m,n]]
R3==Simplify[PB[m -1,n]-PR[m,n]]

B1==Simplify[PR[m+1,n+1]-PB[m,n]]
B2==Simplify[PR[m ,n+1]-PB[m,n]]
B3==Simplify[PR[m +1,n]-PB[m,n]]

Here are diagramatic representations of those results:

({{•, •, •, •}, {•, R1, R3, •}, {•, R2, PR[m,n], •}, {•, •, •, •}})

({{•, •, •, •}, {•, PB[m,n], B2, •}, {•, B3, B1, •}, {•, •, •, •}})

Placing hexagonal node coordinates in linear array

This is the critical first step toward construction of the stochastic vectors and Markov matrices that enter into description of the averaged properties of random walks on hexagonal graphs.

To avoid the extraneous complications that arise when the {m, n}-indices are allowed to assume negative values I will assume that {m, n} ranges on the upper right quadrant of the unit lattice. Suppose, for example, that m and n range on {1, 2, 3, ... , 9}. We confront then the following hexagonal array of 162 nodes:

RedPopulation2=Flatten[Table[PR[m,n],{m,1,9},{n,1,9}],1];
RedPopFig2=Graphics[{Red, PointSize[0.02],Point[RedPopulation2]}, Axes->True, AxesLabel->{x,y}, Ticks->False, AxesOrigin->{0,0}];
BluePopulation2=Flatten[Table[PB[m,n],{m,1,9},{n,1,9}],1];
BluePopFig2=Graphics[{Blue, PointSize[0.02],Point[BluePopulation2]}, Axes->True, AxesLabel->{x,y}, Ticks->False, AxesOrigin->{0,0}];
Show[{RedPopFig2,BluePopFig2}]

Length[RedPopulation2]
Length[BluePopulation2]

To display the elements of the d × d matrix

as a (d^2)-dimensional column vector I will (since I can think of no more sophisticated/natural procedure) simply stack column upon column, in dictionary order:

Generally

where

𝒥[m_,n_]:=d(m-1)+n

as I verify:

𝒥[1,1]//Expand
𝒥[1,2]//Expand
𝒥[1,d]//Expand
𝒥[2,1]//Expand
𝒥[2,2]//Expand
𝒥[2,d]//Expand
𝒥[3,1]//Expand
𝒥[3,2]//Expand
𝒥[3,d]//Expand
𝒥[d,1]//Expand
𝒥[d,2]//Expand
𝒥[d,d]//Expand

Table[(Subscript[a, m,n]),{m,1,9},{n,1,9}]//MatrixForm

To demonstrate (in this illustrative context) the success of

I assign to d the following value

d=9;

and obtain

Flatten[Table[(Subscript[e, 𝒥[m,n]]),{m,1,9},{n,1,9}],1]

To proceed in the reverse direction

(Subscript[e, 𝒥]) ⟶ (Subscript[a, m[𝒥],n[𝒥]])

I use these number-theoretic commands

?Floor

?Mod

to (after a little experimentation) construct

m[x_]:=Floor[((x-1)/d)]+1
n[x_]:=Which[Mod[x,d]≠0,Mod[x,d],Mod[x,d]==0,9]

Here follow series of tests

Table[m[k],{k,1,9}]
Table[n[k],{k,1,9}]

Table[m[k],{k,10,18}]
Table[n[k],{k,10,18}]

Table[m[k],{k,5d+1,6d}]
Table[n[k],{k,5d+1,6d}]

𝒥[RandomInteger[{1,9}],RandomInteger[{1,9}]]

m[48]
n[48]

𝒥[6,3]

culminating in this grand demonstration:

Table[(Subscript[a, m[j],n[j]]),{j,1,(d^2)}]

The hexagonal array

presents a duplex set of doubly-indexed nodal points:

When those (d^2)+(d^2) elements are deployed in columnar array we end up (in the case d = 9) with 162-dimensional column vectors that are to be multiplied by 162 ×162 Markov matrices that have (even though d = 9 constructs a very small sheet of graphene) already this very large number of elements:

(((2×(9^2)))^2)

I elect to deploy those 2(d^2) in a columnar array of this design:

We now have

where as before

It becomes natural at this point to modify that definition:

𝒥a[m_,n_]:=d(m-1)+n
𝒥b[m_,n_]:=(d^2)+d(m-1)+n

𝒥a[1,1]
𝒥a[9,9]
𝒥b[1,1]
𝒥b[9,9]

m[x_]:=Floor[((x-1)/d)]+1
n[x_]:=Which[Mod[x,d]≠0,Mod[x,d],Mod[x,d]==0,9]

e[j_]:=Which[j≤(d^2),(Subscript[a, m[j],n[j]]),j>(d^2),(Subscript[b, m[j-(d^2)],n[j-(d^2)]])]

Table[e[k],{k,1,18}]

Table[e[k],{k,72,90}]

Table[e[k],{k,145,162}]

Structure of the adjacency matrix

Since stochastic vectors have the structure

and on a honeycomb no node is connected to a node of the same color, the adjacency matrix will have the form

where 𝕆 is the (d^2)×(d^2) zero matrix and 𝕊 is a sparce (d^2)×(d^2) matrix with 3 (but occasionally only 2) 1s on each row and all other elements 0. Which means that about

of the elements in 𝕊 are 0s. In the case d = 9 this means that about 96% of the elements in 𝕊 are 0s. The following table shows how that % increases as the value of d ascends:

Table[{d,N[100 (((d^2)-3)/(d^2))]},{d,10,100,10}]//TableForm

We have discovered that the blue nearest neighbors of the red node

(Subscript[{m,n}, red])

lie at

(Subscript[{m-1,n-1}, blue]), (Subscript[{m,n-1}, blue]), (Subscript[{m-1,n}, blue])

we want to express this in terms of the 𝒥 coordinates that range on {1,... ,2(d^2)}.

Clear[d]
𝒥a[m,n]
𝒥b[m,n]

From

𝒥b[m-1,n-1]-𝒥a[m,n]//Expand
𝒥b[m,n-1]-𝒥a[m,n]//Expand
𝒥b[m-1,n]-𝒥a[m,n]//Expand

we obtain

𝒥a[m,n]-1-d+(d^2)==𝒥b[m-1,n-1]//Simplify
𝒥a[m,n]-1+(d^2)==𝒥b[m,n-1]//Simplify
𝒥a[m,n]-d+(d^2)==𝒥b[m-1,n]//Simplify

As m and n range on {1, 2, ... , 9} we see that

run though the expected values

d=9;

Table[𝒥a[m,n],{m,1,9},{n,1,9}]//MatrixForm

Table[𝒥b[m,n],{m,1,9},{n,1,9}]//MatrixForm

but that some of 𝒥-coordinates of the expected blue nearest neighbors of red nodes

lie outside the accepted range {82 ≤ 𝒥 ≤ 162}, as I demonstrate by means of a single example:

Table[𝒥b[m-1,n-1],{m,1,9},{n,1,9}]//MatrixForm

I find by experimentation that to expunge those terms we must arrange for

never to assume negative values. That done, we get

Table[𝒥b[m-1,n-1],{m,2,9},{n,2,9}]//MatrixForm

Table[𝒥b[m,n-1],{m,1,9},{n,2,9}]//MatrixForm

Table[𝒥b[m-1,n],{m,2,9},{n,1,9}]//MatrixForm

To accomplish that snipping in 𝒥-language I introduce a snip factor constructed from

?UnitStep

I proceed from

Clear[d]

𝒥b[m-1,n-1]==𝒥a[m,n]-1-d+(d^2)//Simplify
𝒥b[m,n-1]==𝒥a[m,n]-1+(d^2)//Simplify
𝒥b[m-1,n]==𝒥a[m,n]-d+(d^2)//Simplify

to the definitions

(Subscript[J, 1])[x_]:=x-1-d+(d^2)
(Subscript[J, 2])[x_]:=x-1+(d^2)
(Subscript[J, 3])[x_]:=x-d+(d^2)

Introducing the abbreviation

δ[x_,y_]:=KroneckerDelta[x,y]

we reinstate our illustrative d-value

d=9;

and to construct the (𝒥^th) row of the adjacency matrix 𝔸 command

The (13^th) row, for example, becomes

𝒥=13;Table[UnitStep[j-82](δ[j,(Subscript[J, 1])[𝒥]]+δ[j,(Subscript[J, 2])[𝒥]]+δ[j,(Subscript[J, 3])[𝒥]]),{j,1,162}]

in which

Total[%]

elements are 1, the other 159 are 0. Looking to the 81 × 81-dimensional 𝕊-matrix (within which all the 1s reside), we find the number of 1s in successive rows to be given by

Table[Total[Table[UnitStep[j-82](δ[j,(Subscript[J, 1])[𝒥]]+δ[j,(Subscript[J, 2])[𝒥]]+δ[j,(Subscript[J, 3])[𝒥]]),{j,82,162}]],{𝒥,1,81}]

Length[%]

𝕊=Table[Table[UnitStep[j-82](δ[j,(Subscript[J, 1])[𝒥]]+δ[j,(Subscript[J, 2])[𝒥]]+δ[j,(Subscript[J, 3])[𝒥]]),{j,82,162}],{𝒥,1,81}];

Table[Total[𝕊[[𝒥]]],{𝒥,1,81}]

Det[𝕊]

The idea now is to construct

𝔸=(({{𝕆, 𝕊}, {𝕋, 𝕆}}))

where 𝕆 is the 81-dimensional zero matrix and 𝕋 denotes the transpose of 𝕊. To that end, I construct

ℤ=Table[0,{k,1,81}];
𝕆=DiagonalMatrix[ℤ];

𝕋=Transpose[𝕊];

Dimensions[𝕆]
Dimensions[𝕊]
Dimensions[𝕋]

and borrow from "Linear Algebra of Supersymmetric Matrices" (4-11 August 2011)the ArrayFlatten command, which works this way:

ArrayFlatten[{{(({{0, 0}, {0, 0}})),(({{1, 2}, {3, 4}}))},{(({{5, 6}, {7, 8}})),(({{0, 0}, {0, 0}}))}}]//MatrixForm

We arrive thus at

𝔸=ArrayFlatten[{{𝕆,𝕊},{𝕋,𝕆}}];

Dimensions[𝔸]

But this—by

Det[𝔸]

—cannot be correct. The problem is seen by

Table[Total[𝔸[[k]]],{k,1,162}]

to trace to the fact that the first and last row/column of 𝔸 are all 0, meaning that there is a red node that connects to no blue node, and a blue node that connects to no red node. Which is, in fact, obvious from a figure that has already been displayed twice

and which, as has already beem remarked, refers to a disconnected graph (the top ● and bottom ● are disconnected from the rest of the graph).

To address this problem I propose tentatively to
   • reduce the dimension of 𝕆  from 81 to 80;
   • delete the first row and last column from 𝕊 (so that it too will become 80 × 80).

𝕫=Table[0,{k,1,80}];
𝕠=DiagonalMatrix[𝕫];
Dimensions[𝕠]

𝕤=Table[Table[UnitStep[j-81](δ[j,(Subscript[J, 1])[𝒥]]+δ[j,(Subscript[J, 2])[𝒥]]+δ[j,(Subscript[J, 3])[𝒥]]),{j,81,160}],{𝒥,2,81}];
Dimensions[𝕤]
Table[Total[𝕤[[𝒥]]],{𝒥,1,80}]
Det[𝕤]

𝕥=Transpose[𝕤];

𝕒=ArrayFlatten[{{𝕠,𝕤},{𝕥,𝕠}}];
Dimensions[𝕒]
Det[𝕒]
Table[Total[𝕒[[k]]],{k,1,160}]

REMARK: The palindromic structure of that list—in which 1s, 2s and 3s occur in the following numbers

Tally[%]

—is reminiscent of the design of the precediong figure.

Construction and properties of the associated Markov matrix

((𝕒[[1]])/(Total[𝕒[[1]]]))

((𝕒[[8]])/(Total[𝕒[[8]]]))

((𝕒[[10]])/(Total[𝕒[[10]]]))

𝕄=Table[((𝕒[[k]])/(Total[𝕒[[k]]])),{k,1,160}];

Dimensions[𝕄]

Det[𝕄]

N[%]

The following calculation takes maybe two minutes to run

Eigenvalues[𝕄][[1]]

and produces a mysterious minus sign which presumably has to do with the fact that the honeycomb graph is a "blinker." In support of that suspicion, we compute

Eigenvalues[𝕄][[2]]

I entered the following command shortly before heading off to attend Max Schlosshauer's seminar (3:40 PM on 25 January 2012). By 10:00 AM the next morning Mathematica still had not produced a result.

Eigenvectors[𝕄][[1]]

Mathematica is, however, able to evaluate moderately large powers of 𝕄 in a few seconds. Here I look to the top rows of (𝕄^20) and (𝕄^21):

MatrixPower[𝕄,20][[1]]//N

MatrixPower[𝕄,21][[1]]//N

Those results conform to the general propositions that

and place us in position to guess the structure of the eigenvectors that associate with the eigenvalues ±1. From the following "red/blue lists"

r=Join[Table[1,{j,1,80}],Table[0,{k,1,80}]]
b=Join[Table[0,{j,1,80}],Table[1,{k,1,80}]]

construct column vectors

𝕣=Transpose[{r}];
𝕓=Transpose[{b}];

and observe that

𝕄.𝕣==𝕓
𝕄.𝕓==𝕣

𝕄.(𝕣+𝕓)==+1(𝕣+𝕓)

𝕄.(𝕣-𝕓)==-1(𝕣-𝕓)

These results speak clearly of "blinking equilibration."

Concluding remarks

The "corrugated graphene problem" would lead to Markov matrices derived from weighted adjacency matrices. On finite honeycombs (however large) we expect to be led asymptotically to a distorted version of blilnking equilibration. But we cannot expect in such cases to guess the structure of the relevant eigenvectors (those associated with eigenvalues ±1). And to obtain them by computation would—if the lattice is of physically interesting size—require a supercomputer.

I am reminded that the computational requirements of lattice gauge theory

Hyperlink["Lattice Gauge Theory","http://en.wikipedia.org/wiki/Lattice_gauge_theory"]

served historically as a major stimulus toward the development of supercomputers.

I have not attempted to introduce site-specific weighting, whether along lines suggested by naive corrugation models or otherwise. I think I see how to do so (by adapting techniques explored in an earlier notebook), but in view of the computational problem encountered in the preceding section I see no reason to introduce that added degree of complexity, since I would be computationally prevented from exploring its consequences.